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This  report  is  divided  into  four  parts:  Part-I  provides  an  introduction  to 
plunging  jet  phenomena  and  summarizes  the  major  findings  of  our  previous 
research.  Part-II  presents  a  summary  of  the  major  experimental  findings 
obtained  during  at  Rensselaer.  Part-Ill  describes  two  proposed  mechanisms 
responsibles  for  air  entrainment.  Part-IV  presents  the  results  of  two  different 
numerical  calculations,  one  using  a  Euler-Lagrangian  approach  and  the  other 
based  on  solving  a  full  two-fluid  model  with  a  CFD  code.  Finally  Part-V  gives 
conclusions  and  recommendations  for  future  research. 
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PART-I:  INTRODUCTION 

The  gas  entrained  by  a  plunging  liquid  jet  and  the  resultant  two-phase  jet 
dispersion  occur  in  many  problems  of  practical  interest.  In  particular,  the  air 
entrainment  process  due  to  the  breaking  bow  waves  of  surface  ships  may  cause 
long  (ie,  up  to  5  km  in  length)  wakes.  Naturally  easily  detectable  wakes  are 
undesirable  for  naval  warships.  The  ecological  balance  in  lakes  and  oceans  is 
also  dependent  on  the  amount  of  dissolved  oxygen.  A  good  understanding  of  the 
air  carryunder  and  bubble  dispersion  process  associated  with  a  plunging  liquid 
jet  is  vital  if  one  is  to  be  able  to  quantify  such  diverse  phenomena  as  sea  surface 
chemistry,  the  meteorological  significance  of  (breaking)  ocean  waves,  the 
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performance  of  certain  type  of  chemical  reactors,  the  "greenhouse"  effect  (ie,  the 
absorption  of  CO2  by  the  oceans),  and  a  number  of  other  important  maritime- 
related  applications.  Significantly  the  absorption  of  greenhouse  gases  into  the 
ocean  has  been  hypothesized  to  be  highly  dependent  upon  the  air  carryunder  that 
occurs  due  to  breaking  waves.  This  process  can  be  approximated  with  a  plunging 
liquid  jet  (Monahan  [1991],  Kerman  [1984]).  In  addition,  the  air  carryunder  that 
occurs  at  most  hydraulic  structures  in  rivers  is  primarily  responsible  for  the 
large  air/water  mass  transfer  that  is  associated  with  these  structures  (Avery  and 
Novak  [1978]). 


Also,  air  entrainment  plays  an  important  role  in  slug  flow  phenomena  in 
conduits.  Indeed  the  liquid  film  surrounding  the  Taylor  bubble  has  a  mean  flow 
in  the  opposite  direction  from  the  Taylor  bubble.  This  liquid  forms  a  type  of 
plunging  jet  that  produces  a  surface  depression  in  the  rear  part  of  the  Taylor 
bubble.  When  the  annular  liquid  jet  exceeds  a  critical  velocity,  the  plunging  liquid 
jet  entrains  into  the  liquid  slugs  small  bubbles  from  the  air  in  the  Taylor  bubble. 
These  bubbles  follow  and  may  coalesce  with  the  Taylor  bubbles  . 

A  number  of  prior  studies  have  been  performed  in  which  axisymmetric 
plunging  liquid  jets  have  been  used  to  investigate  the  air  carryunder  process. 
These  include  the  work  of  Lin  &  Donnelly  [1966],  Burgess  et  al  [1972],  Van  De 
Sande  &  Smith  [1973],  Koga  [1982],  McKeogh  &  Ervine  [1981],  and  Detsch  & 
Sharma  [1991],  Ohkawa  et  al  [1986]  Ervine  et  al  [1980],  McKeogh  &  Elsaway  [1980], 
Ervine  &  Falvey  [1987],  Blanchard  &  Cipriano  [1981],  and  Sene  [1988]. 

The  primary  objective  of  the  experimental  research  was  to  obtain  detailed 
local  data  in  the  two-phase  flow  region  of  a  plunging  liquid  jet  that  could  be  used 
to  validate  multidimensional  two  fluid  model  CFD  calculations.  We  have 
measured  the  local  (turbulent)  velocity  of  the  liquid,  the  velocity  of  the  gas  bubbles, 
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the  bubble  size  distribution  for  bubbles  smaller  than  about  1.0  mm  and  the  void 
fraction  of  the  gas  phase.  The  combined  probability  density  function  of  the  bubble 
size  and  velocity  has  not  been  measured  before.  The  probability  density  function 
for  bubble  diameter  is  needed  to  compute  the  interfacial  area  density,  an 
important  parameter  which  helps  determine  the  mass,  momentum  and  energy 
transfer  characteristics  of  a  two-phase  flow.  The  probability  density  function  of 
the  bubble  velocity  is  needed  to  estimate  the  total  time  that  the  bubbles  remain 
submerged  and  therefore  able  to  transfer  mass,  momentum  and  energy. 

The  spreading  of  a  two-phase  jet  involves  the  interaction  between  the  liquid 
turbulence  and  the  bubbles.  This  problem  can  not  be  solved  analytically.  However 
surface  depression  and  air  entrainment  can  be  modeled  using  simplifications, 
such  as  inviscid-irrotational  flow  theory. 

Recently  Lezzi  and  Prosperetti  (1991)  have  proposed  that  the  instability 
responsible  for  the  air  entrainment  was  caused  by  the  gas  viscosity.  In 
particular,  they  studied  the  linear  stability  of  a  vertical  film  of  a  viscous  gas 
bounded  by  an  inviscid  liquid  in  uniform  motion  on  one  side,  and  by  inviscid  liquid 
at  rest  on  the  other  side.  In  the  long  wavelength  asymptotic  limit,  the  gas  film 
behaves  as  a  single  surface  and  thus  they  obtained  the  classical  Helmholtz  result 
when  only  a  surface  is  present.  They  also  obtained  the  marginal  stability 
boundary  with  the  gas  gap  width,  8,  as  the  control  parameter  (i.e.,  they 
numerically  compute  for  a  given  8  the  range  of  wavelength  that  makes  the  system 
unstable). 

For  the  computation  of  jet  dispersion  it  is  important  to  appropriately  model 
the  turbulent  intensity  of  the  continuous  phase.  Single-phase  turbulent  jets 
represent  an  important  class  of  free  shear  flows  that  have  been  studied  in  the  past 
to  develop  and  test  turbulence  models  [Abramovich,  1963].  More  recently, 
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turbulent  jets  have  been  evaluated  numerically  using  computational  fluid 
dynamic  (CFD)  techniques  for  various  turbulence  models.  Rodi  [1984]  presented 
results  using  the  classical  k-e  model  of  Gibson  &  Launder  [1976],  and  showed  that 
k-e  models  may  not  accurately  predict  jet  spreading.  Rodi  [1984]  proposed  that  the 
constant  in  the  model  for  turbulent  viscosity  was  really  a  function  of  the  ratio 

between  turbulent  production  and  dissipation.  This  involved  the  development  of  a 
new  function  which  produced  better  results.  Sini  and  Dekeyser  [1987]  solved  the 
single-phase  turbulent  jet  using  Rodi’s  k-e  model  [Rodi,  1984].  This  model 
compared  favorably  with  the  experimental  results  of  the  single-phase  turbulent 
jet  as  well  as  with  other  more  detailed  algebraic  stress  models.  Hence  it  appears 
that  in  some  cases  turbulent  nonisotropy  is  not  important  and  need  not  be 
modeled. 

Significantly,  it  has  been  found  that  single-phase  turbulent  jet  data  can  be 
used  for  the  assessment  of  turbulence  models  because  one  does  not  have  to 
constitute  complicated  turbulent  closure  laws  near  solid  (no  slip)  boundaries. 
Indeed,  due  to  the  absence  of  walls  and  the  associated  shear  boundary  conditions 
the  turbulent  jet  is  probably  the  simplest  non-trivial  case  to  analyze. 
Interestingly,  the  same  conclusions  can  be  reached  for  a  two-phase  turbulent  jet. 

In  most  of  the  previously  mentioned  research,  the  flow  field  was  considered 
to  be  thin  in  the  lateral  direction  and  the  flows  were  characterized  by  a  relatively 
small  lateral  velocity  when  compared  to  the  streamwise  velocity.  This  is 
equivalent  to  considering  the  flow  (ie,  the  liquid  jet)  to  be  a  boundary  layer.  Hence, 
the  flow  is  often  assumed  to  have  a  negligible  pressure  gradient  in  the  lateral 
direction,  and  thus  the  pressure  in  the  boundary  layer  can  be  imposed  by  the 
external  flow.  From  the  mathematical  point  of  view,  this  approximation  changes 
the  problem  from  an  elliptic  to  parabolic  one  in  the  streamwise  direction.  From 
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the  computational  point  of  view,  when  using  a  boundary  layer  approximation  the 
pressure  distribution  has  to  be  prescribed.  That  is,  the  static  pressure  as  a 
function  of  the  streamwise  coordinate  has  to  be  known  and  supplied  to  the 
computational  fluid  dynamic  (CFD)  code. 

Depending  on  the  variables  of  interest,  approximating  the  jet  as  a  boundary 
layer  may  be  a  reasonable  approximation.  However,  the  lateral  velocity  at  the 
edge  of  the  jet  (i.e.,  the  entrainment  velocity)  may  be  much  larger  than  the 
streamwise  velocity.  Moreover,  for  a  planar  jet  the  entrainment  velocity  remains 
finite  as  the  integration  domain  in  the  lateral  direction  is  enlarged. 

In  this  work  we  did  not  make  the  boundary  layer  approximation.  Rather, 
we  assumed  an  elliptic  problem  for  both  the  gas  and  the  liquid  phases.  Solving 
the  partial  differential  equations  as  an  elliptic  system  increases  the  complexity  of 
the  numerical  problem  but  provides  more  detailed  and  accurate  information  on 
the  flow  field  than  a  parabolic  scheme.  One  of  the  advantages  of  the  elliptic 
solution  is  that  the  streamwise  pressure  distribution  does  not  have  to  be  provided. 
This  has  two  effects,  first  we  are  able  to  compute  recirculation,  and  second,  we 
are  able  to  calculate  the  buoyancy-induced  reversal  of  the  bubble  motion.  These 
two  important  effects  cannot  be  computed  using  a  parabolic  approach. 

In  parabolic  single-phase  jet  calculations  using  a  k-e  model,  Sini  & 
Dekeyser  [1987],  and  Hossain  &  Rodi  [1982]  found  satisfactory  agreement  between 
calculations  and  experiments,  except  for  a  small  region  near  the  jet  inlet.  Their 
good  results  are  due  in  part  to  the  fact  that  they  analyzed  a  free  jet  which  was  only 
weakly  nonisotropic.  The  next  level  of  complexity  for  simulating  turbulent  two- 
phase  flows  is  the  use  of  an  Algebraic  Stress  Model  (ASM).  Several  different 
models  have  been  proposed  in  the  literature.  For  the  particular  case  of  a  planar 
jet,  performance  of  the  ASM  is  similar  to  the  k-e  model.  One  might  expect  that 
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the  nonisotropic  ASM  model  would  produce  better  results  than  the  isotropic  k-e 
model  for  the  pressure  distribution.  However,  when  we  used  the  ASM  proposed 
by  Gibson  &  Launder  [1976]  for  the  evaluation  of  a  planar  jet  the  computed 
streamwise-pressure  distribution  [Bergstrom,  1992]  was  virtually  the  same  as 
when  a  k-e  model  was  used. 

Thus,  for  simplicity,  we  have  used  a  k-e  model  in  the  computations 
presented  in  this  report.  We  note  that  the  turbulence  present  in  the  liquid  has  two 
components  in  a  two-phase  jet.  One  component,  the  shear-induced  turbulence,  is 
due  to  viscosity  and  it  is  present  in  both  single  and  two-phase  flows.  The  other 
component  is  the  bubble-induced  turbulence  due  to  slip  between  the  bubbles  and 
the  surrounding  liquid,  and  it  only  occurs  in  two-phase  flows. 

PART-II  EXPERIMENTAL  RESULTS 

As  shown  in  Figure- 1,  a  converging  nozzle  oriented  vertically  produced  an 
axisymmetric  liquid  (ie,  water)  jet.  This  jet  impacted  at  90°  a  pool  of  water  and, 
when  a  threshold  velocity  was  exceeded,  it  was  observed  that  the  plunging  liquid 
jet  caused  air  entrainment.  In  agreement  with  the  observations  of  McKeogh  & 
Ervine  [1981],  different  two-phase  jet  characteristics  were  noted,  depending  on  the 
turbulence  intensity  of  the  plunging  liquid  jet.  For  a  laminar  liquid  jet  (ie,  one 
having  a  turbulence  intensity  less  than  about  0.8%)  the  diameter  of  the  entrained 
bubbles  were  in  the  range  15-300  |im.  On  the  other  hand,  for  liquid  jet  turbulence 
intensity  of  about  3%,  the  entrained  bubbles  had  diameters  in  the  range  of  1-3  mm. 
Our  definition  of  a  rough  and  smooth  jet  should  be  considered  only  as  an  easy  way 
to  refer  to  one  particular  turbulence  intensity.  In  particular,  rough  (smooth) 
means  a  turbulence  intensity  of  about  3%  (.8%).  Another  candidate  for  the 
distinction  was  the  Reynolds  number.  However,  this  is  not  appropriate  because 
after  the  contraction,  with  the  same  Reynolds  number  we  may  have  different 
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turbulence  intensity  and  these  different  turbulence  intensities  generate 
qualitatively  different  two-phase  flows. 

Figure- 1  also  shows  a  schematic  of  the  test  loop.  A  screw  pump  was  used  to 
force  the  water  through  the  nozzle  as  well  through  a  bypass.  The  pump  had  a 
speed  controller  which  was  used  to  make  the  coarse  control  of  the  liquid  flow  rate 
through  the  nozzle.  In  the  bypass  a  valve  was  used  for  the  fine  control  of  liquid 
flow  rate.  In  order  to  damp  out  any  flow  oscillations,  an  accumulator  was  placed 
on  the  discharge  side  of  the  pump. 

The  acrylic  conical  nozzle,  shown  schematically  in  Figure-2,  consisted  of 
an  arrangement  of  honeycombs  and  screens  followed  by  a  smooth  contraction.  In 
this  way  the  turbulence  level  of  the  liquid  jet  could  be  parametrically  varied.  The 
exit  diameter  of  the  nozzle  was  5.1  mm,  and  this  produced  a  liquid  jet  of  about  the 
same  diameter.  The  acrylic  tank  which  contained  the  water  pool  had 
dimensions,  0.914  x  0.916  x  1.465  =  1.265  m3.  The  suction  of  the  tank  was  put  as 
far  from  the  liquid  jet  impact  point  as  possible  in  order  to  minimize  the  influence 
of  this  flow  on  the  two-phase  jet. 

A  DANTEC  Fiber-Flow  Laser  Doppler  Anemometer  (LDA)  system  was  used 
to  nonintrusively  measure  the  liquid  and  gas  velocities  (both  the  mean  and 
fluctuations).  This  system  consisted  of  submersible  transmitting  and  receiving 
optics. 

The  receiving  optics  used  for  the  smooth  jet  employed  a  Fiber  Phase  Doppler 
Anemometer  (FPDA)  system  with  600  mm  focal  length  lenses  and  a  special 
aperture  plate  to  maximize  the  bubble  size  range.  The  axial  velocity  of  the  liquid 
jet  was  used  as  the  master  signal  for  data  collection.  The  collected  light  was 
transmitted  through  three  optical  fibers  to  a  special  FPDA  device  having  three 
photodetectors  for  bubble  size  measurement. 
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The  signals  collected  by  the  AT  computer  consisted  of  the: 

-  arrival  time  of  the  particles 

-  transit  time  of  the  particles 

-  velocity  of  the  particles 

-  equivalent  diameter  of  the  particles 

Two  different  methods  were  used  for  the  measurei  it  of  void  fraction  in 
the  two-phase  jet:  a  KfK  impedance  probe  and  the  particle  time  fraction  from  the 
FPDA.  The  impedance  probe  consists  of  two  electrically  isolated  electrodes;  one  at 
the  tip  of  the  probe  and  another  downstream  electrode  which  was  always  in 
contact  with  the  liquid  in  the  pool.  The  liquid  (ie,  water)  had  a  relatively  high 
electrical  conductivity  and  thus  when  the  tip  was  in  contact  with  the  liquid  a 
relatively  high  current  flows  through  to  a  Wheatstone  bridge  circuit.  The 
difference  between  the  conductivities  of  the  liquid  and  gas  phases  produced  a 
different  signal  depending  on  whether  there  is  liquid  or  gas  present  at  the  tip  of 
the  probe.  The  active  element  of  the  probe’s  tip  was  150  pm  in  diameter  and  it  was 
calibrated  with  bubbles  having  diameters  in  the  1-3  mm  range.  This  type  of  probe 
is  a  standard  tool  used  for  measuring  local  void  fraction  in  air/water  bubbly  flows 
[Hewitt,  1978]. 

The  KfK  impedance  probe  was  used  to  measure  the  void  fraction  in  the 
rough  jet  (which  had  bubbles  of  diameter  in  the  range  1-3  mm)  because  the  size  of 
bubbles  produced  was  out  of  range  of  the  FPDA  (ie,  the  air  bubbles  were  too  large 
for  the  lens  size  used).  In  contrast,  for  the  measurement  of  void  fraction  when  a 
smooth  liquid  jet  was  tested,  the  impedance  probe  could  not  be  used  because  the 
size  of  the  bubble  were  the  same  order  of  magnitude  as  the  size  of  the  tip  (ie,  the 
air  bubbles  were  quite  small).  However,  the  FPDA  could  be,  and  was,  used  to 
measure  the  size  distribution  of  the  bubbles  in  this  case. 


8 


The  FPDA  was  calibrated  using  a  suspension  of  polystyrene  particles 
which  had  a  diameter  of  9.5  pm  ±  0.5  pm,  and  a  rotating  steel  ball  of  diameter 
0.4  mm. 

The  LDA/FPDA  system  and  the  KfK  impedance  probe  were  mounted  on  a 
Benjamin  Systems  three-dimensional  traversing  mechanism  having  a  1  pm 
positioning  resolution.  The  tip  of  the  KfK  impedance  probe  was  positioned  0.3  mm 
under  the  measurement  volume  of  the  LDA/FPDA  system  for  void  fraction 
measurements  when  a  rough  jet  liquid  was  tested. 

The  turbulence  intensity  of  the  liquid  jet  at  the  nozzle  exit  was  found  to  be 
one  of  the  most  important  parameters  affecting  jet  roughness  and  the  size  of  the 
bubbles  entrained  by  the  plunging  liquid  jet.  An  arrangement  of  honeycombs  and 
screens  were  used  to  control  the  turbulence  intensity  of  the  flow  entering  the 
conical  nozzle. 

Figure-3  depicts  a  contour  plot  of  the  two-dimensional  probability  density 
function  of  the  particle  diameters  (dp)  and  axial  velocities  (uz).  Quantitatively  the 
most  probable  value  of  peak  #1  was  at,  dp  =  5  pm,  uz  =  4.05  m/s,  and  the  most 
probable  value  of  peak  #2  was  at  dp  =  125  pm,  uz  =  3.5  m/s. 

Figure-3  indicates  that  the  velocity  of  the  bubble  was  not  dependent  on  its 
size.  If  the  velocity  of  the  bubbles  had  changed  with  size  we  would  see  the  iso¬ 
count  curves  with  their  principal  axes  forming  an  angle  with  the  horizontal. 
There  is  no  such  trend  seen  in  Figure-3. 

Figure-4  shows  the  liquid  velocity  histogram  at  the  centerline  of  the  liquid 
jet  for  a  flow  rate  of  w  =  0.144  kg/s,  a  distance  from  the  nozzle  to  the  undisturbed 
pool  surface  of  h  =  30.0  ±  0.3  mm,  and  a  distance  from  the  pool  surface  to  the 
measurement  volume  of  z  =  33  mm.  The  mean  axial  velocity  is,  ui  =  4.96  m/s. 
One  of  the  main  differences  between  a  rough  liquid  jet  and  a  smooth  jet  is  that  in 
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the  latter  case  the  liquid  flow  field  is  practically  unaffected  by  the  bubbles  while  in 
the  former,  the  bubbles  are  much  larger,  thus  the  discrete  phase  increases  the 
continuous  phase's  turbulence  intensity.  This  also  increases  phasic  momentum 
exchange  resulting  in  greater  dispersion  of  the  two-phase  jet  and  a  lower  velocity 
of  the  liquid  velocity. 

Figure-5  presents  the  liquid  and  gas  velocity  as  a  function  of  radial 
position  (r)  for  h  =  17.3  mm,  wj  =  0.125  kg/s  and  z  =  50.0  mm.  We  see  that  at  the 
edge  of  the  spreading  two  phase  jet  that  the  gas  (bubble)  velocity  is  negative, 
indicating  buoyancy-driven  countercurrent  flow.  It  was  also  found  that  the  two- 
phase  jet  was  more  dispersed  than  the  corresponding  single-phase  flow  case  and 
that  the  turbulence  intensity  was  higher.  The  turbulence  enhancement  is  due  to 
bubble-induced  turbulence.  In  this  case  the  bubble-induced  turbulence  accounts 
for  about  30%  of  the  total  turbulence  level. 

As  noted  previously,  when  the  liquid  jet  impacts  the  pool  surface,  air 
entrainment  occurs  around  the  jet's  circumference.  In  Figure-6a  the  measured 
local  void  fraction  is  presented  as  a  function  of  r  for  z=  1  mm  (i.e.  with  the  probe 
1  mm  under  the  undisturbed  liquid  level).  We  see  that  the  void  fraction  has  a 
maximum  at  r  =  djet,/2  =  2.5  mm.  Obviously,  the  air  entrainment  process  is 
responsible  for  this  effect.  In  the  high  speed  video  visualization  of  these 
experiments  it  was  rare  to  observe  bubbles  at  the  liquid  jet’s  centerline  for  z  <  10 
mm.  However,  once  the  air  was  entrained,  dispersion  of  the  gas  phase  occurred 

as  z  was  increased.  Figure-6b  shows  how  the  void  peaks  in  Figure-6a  was 
dispersed  at  z  =  18  mm.  We  see  that  the  maximum  now  occurs  at  r  =  5  mm. 

Moreover,  we  see  that  there  is  significant  void  fraction  at  r  =  0  (ie,  the  jet’s 
centerline)  because  of  the  void  dispersion  process.  Figure-6c  shows  the  void 
fraction  profile  at  z  =  43  mm.  Significantly,  the  curve  now  has  a  maximum  at  the 
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centerline  of  the  jet  (r  =  0).  Again,  this  is  a  direct  result  of  the  void  dispersion 
process  in  the  two-phase  jet. 

Figures-7a,  7b  and  7c  present  the  axial  liquid  velocity  (uj)  as  a  function  of 
the  radial  distance  r  for  three  different  distances  (z)  from  the  undisturbed  pool 
level.  Figure-7  a  is  a  measurement  of  the  liquid  velocity  profile  at  the  exit  of  the 
nozzle.  This  flat  profile  is  characteristic  of  the  potential  flow  which  exited  from 
the  nozzle. 

Figures-7b  and  7c  show  that  the  liquid  velocity  curves  have  a  maximum  at 
the  centerline  for  all  z.  Comparing  these  two  profiles  one  can  easily  see  how  the 
two-phase  jet  spreads. 

Figure-8  presents  the  liquid  and  gas  phasic  mean  velocities  as  a  function  of 
the  lateral  position  (y)  for  a  distance  between  the  nozzle  exit  and  the  undisturbed 
pool  level  of  h=8.5mm,  a  liquid  jet  flow  rate  of,  w=1.8  kg/s,  and  a  submergence  of, 
z=31mm.  We  see  that  at  the  centerline  the  relative  velocity  is  approximately 
0.28m/s.  This  is  very  close  to  the  terminal  velocity  in  water  of  a  gas  bubble  having 
a  3  mm  diameter.  The  ensemble  averaged  two-fluid  model  presented  in  Part-IV 
predicts  that  at  a  symmetry  plane  (where  all  lateral  velocity  gradients  are  zero 
due  to  symmetry)  the  relative  velocity  in  the  axial  direction  is  determined  by  a 
balance  between  buoyancy  and  drag  (ie,  the  same  forces  that  determines  the 
terminal  velocity).  Our  experimental  findings  support  this  conclusion. 

We  see  in  Figure-8  that  at  the  edge  of  the  spreading  jet  the  gas  velocity  is 
negative,  indicating  buoyancy-driven  countercurrent  flow.  It  was  also  found  that 
the  two-phase  jet  was  more  dispersed  than  the  corresponding  single-phase  flow 
case.  This  is  presumably  due  to  the  fact  that  the  presence  of  the  gas  phase 
obstructs  a  fraction  of  the  flow  area  for  the  liquid  producing  an  additional 
acceleration  of  the  liquid  and  therefore  a  larger  momentum  interchange  with  the 


liquid  pool.  We  also  observe  that  the  turbulence  instensity  was  higher  than  in  the 
single-phase  flow  due  to  bubble-induced  pseudoturbulence. 

Figure-9  shows  the  void  fraction  as  a  function  of  the  lateral  position  (y)  for 
three  different  axial  location  (z).  The  spreading  of  the  dispersed  phase  (ie,  the 
bubbles)  can  be  easily  seen  in  these  plots.  Figure-9  also  shows  the  curve  for 
z=lmm  (ie,  the  probe  1mm  under  the  surface)  where  one  can  see  a  void  fraction 
peak  at  about  y=3.5  mm.  This  peak  is  close  to  the  liquid  jet's  edge  (y=2.05mm)  and 
most  of  the  air  is  entrained  in  the  peak  neigborhood.  For  y=0mm  (centerplane)  at 
z=lmm  the  void  fraction  has  a  minimum.  This  point  corresponds  to  the  liquid 
jet's  center  impacting  the  probe  and  threrefore  the  bubble  population  there  is 
very  smallsince  little  lateral  dispersion  has  ocurred.  Finally  Figure-9  shows  how 
the  void  fraction  peak  off  the  centerline  is  reduced  as  we  move  down  into  the  pool 
and  also  shows  significant  void  dispersion.  Figure  10  depicts  the  centerline  void 
fraction  as  a  function  of  the  axial  position.  We  can  also  see  here  the  local  void 
fraction  builds  up. 

PART-IH:  ANALYTICAL  RESULTS 

The  objective  of  this  section  is  to  evaluate  the  induced  surface  depression 

caused  by  a  plunging  plane  liquid  jet.  We  have  assumed  that  both  fluids  are 

inviscid  and  irrotational,  and  that  the  gas  region  is  at  constant  pressure.  Hence 

the  appropriate  equation  for  this  problem  is  Laplace's  equation  for  the  liquid  and 

the  associated  interfacial  jump  condition.  The  problem  may  be  described  by  two 

12 

2  P*  xo  P/gx0 

parameters,  the  Weber  number,  We  = - - - ,  and  the  Bond  number,  Bo=  — — — 

,  where  x0  is  the  half- width  of  the  plunging  liquid  jet,  v(  is  the  velocity  of  the  liquid 

jet,  Pf  is  the  density  of  the  liquid,  g  is  the  gravitational  acceleration,  and  c  is  the 

surface  tension. 
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We  have  used  non-singular  perturbation  techniques  to  analyze  the  problem. 
That  is,  we  have  expanded  the  solutions  for  relatively  small  We,  and  substituted 
these  expressions  into  Laplace's  equation  and  the  interfacial  jump  condition. 
Equating  terms  of  the  same  order  in  We,  we  obtained  a  recursive  system  of 
equations,  which  were  solved  numerically  up  to  third  order. 

For  the  assumption  of  irrotational,  inviscid  flow  the  governing  equations  for 
the  liquid  field  are, 

V2y  =  0  (1) 

The  two  velocity  components  are  then  given  as, 

(2) 

(3) 


where  y  is  the  stream  function. 
3y 


u  = 


V  = 


ay 

dy 

3x 


where  u  and  v  are  the  velocity  components  along  the  horizontal  x  and  vertical  y 

axis,  respectively.  If  we  prescribe  a  value  of  y,  or  its  derivative  normal  to  the 

ay  ,  , 

surface,  pr,  for  every  point  on  the  boundary,  we  have  a  well-posed  problem.  As 
dn 

shown  in  Figure-11,  the  plane  x  =  0  is  a  symmetry  plane.  Thus  the  transverse 
lateral  velocity,  u,  must  vanish  for  every  y  atx  =  0.  That  is,  from  Eq.  (2), 


u(x  =  0  ,y)  =  |^r  (x  =  0  ,y)  =  0 


(4) 


on  the  centerline  of  the  plunging  liquid  jet.  We  know  that  the  axial  velocity  must 
be  v(  at  points  where  the  jet  is  impacting.  From  Figure-11  and  Eq.  (3),  we  see 

that: 
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V  (0  <x  <  x0,y  =0)  =  -||-(0  <x  <  x0,y  =  0)  =  v, 


(5) 


For  x  >  x0,  the  free  surface  position,  T]  (x),  must  be  coincident  with  a  streamline. 

Without  loss  of  generality,  we  can  make  the  free  surface  coincident  with  the 
streamline  ijf  =  0.  Then, 

\/(x,f|(x))  =0  (6) 


We  assume  that  the  pressure  in  the  gas  region  is  constant  and  equal  to  zero.  The 
liquid  pressure  right  under  the  surface  is  related  to  the  curvature  of  the  surface 
by: 


d2* 


1  + 


'dn 

dx 


-^  =  -p(S,Tl(x)> 


(7) 


Bernoulli’s  equation  for  the  liquid  pressure  gives 

-p^  (x,T|(x))  =  |p^  [u2(x,“n(x)) +v2(x,fi(x))]  + p^gii(x)  (8) 

We  may  make  Eqs.  (l)-(8)  nondimensional,  using  x0  as  the  length  scale  and,  v(  as 
the  velocity  scale,  respectively.  Thus  we  have: 


V2V  =  0  (9) 

with  boundary  conditions, 

|y  (x  =  0,  y)  =  0  (10a) 

V  (x,  y->«»)  =  0  (10b) 
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V  (x-*»,y)  =  0 

(10c) 

1-4 

II 

o 

II 

►» 

SH  * 

rt>|ru 

0<x<  1 

(10d) 

y(x,Ti  (x))  =  0 

1  <  X 

(lOe) 

where, 

X  =x/x0 

O 

& 

n 

>» 

.  11  =T1/X0 

u  =u l\t 

< 

it 

<> 

< 

A 

V 

and  v=XVJ 

The  velocity  components  may  now  be  computed  from  the  stream  function  using 
Eqs.  (2)  and  (3)  as, 

u(x,y)  =  (x,y)  (11) 

v(x,y)  = (x,y)  (12) 

We  shall  solve  the  problem  for  small  We  using  a  perturbation  analysis. 
First,  we  may  expand  all  the  dependent  variables  in  terms  of  We  obtaining: 
n 

ti(x)  =  Z  We*  'nj(x)  +  0(Wen+1)  (13a) 

i=l 

n 

v(x,y)  =  I  We*  VjCx.y)  +  0(Wen+1)  (13b) 

i=0 

Substituting  eqs.  (13)  in  eqs.  (10)  and  collecting  terms  of  the  same  order  in  We  we 
obtain  the  following  boundary  conditions  for  x  >  1 , 

yo(x,0)  =  0  (14a) 
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y^x.0)  =  -  u0ti1 


(14b) 


1  2 

V2(x,0)  =  -I  U1T11  +  u0ti2  +  217^1 


(14c) 


f  l  dUj  2  1  *^no  3^1 

V3(x,0)  =  -\xxQr\3  +  Ulri2  +  u2ti1  +  g-g^i  +  ^^2  +  J  (14d) 

Notice  that  \|f.  has  homogeneous  boundary  conditions  for  i  >  1,  except  for  x  >  1, 

where  its  value  is  given  by  Eqs.  (14). 


Substituting  eqs.  (13)  in  eq.  (7)  and  collecting  terms  of  the  same  order  in  We  we 
obtain  the  following  boundary  conditions  for  x  >  1  , 


=  Bo  Tl 


(14a) 


d2q,  o  2 

-^r  =  (u0+v0)  +  BoTh 

d\  f  fdu  o  1  f^o  ^ 

TT«  2Uo  '5"T,1  +  U1  +2vo  '^Tll  +  Vl  +  BoTfc 


(14b) 


(14c) 


The  boundary  conditions  for  these  ordinary  differential  equations  are: 


Tlj  (x— >°°)  =  0 
drij 

^(x->~)  =  0 


(15a) 


(15b) 


Equation  (15b)  results  in  the  following  ordinary  differential  equation, 
t^'-Boti^-OIo  +  Vq)  (17) 


for  the  interval  1  <  x  <  D  with  the  boundary  conditions, 


ti1(x=D)  =  0 


(18a) 
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(18b) 


dth 

dx 


(x=D)  =  0 


Equation  (17)  has  an  integrable  singular  point  at  which  x=l.  We  have 

numerically  evaluated  Eq.  (17)  using  a  Runge-Kutta  (RK)  algorithm  with  adaptive 
stepsize.  Figure- 12  shows  as  a  function  of  x. 

Figure-13  shows  the  position  of  the  interface  q  as  a  function  of  the  lateral 
position  x  for  values  of  the  Weber  number  in  the  range  [0.1  to  0.8].  We  now 
analyze  the  behavior  of  the  q(x)  as  we  move  from  the  rightmost  position  in  the  plot 
(x=l.l)  to  lower  values  of  x.  For  We  =  0.1,  q(x)  increases  as  x  approaches  the  jet 
impact  point  x  =  1  without  any  local  minimum  (ie,  no  surface  depression  is 

present).  For  We  =  0.8,  on  the  other  hand,  as  x  is  decreased  we  first  have  a  local 
maximum  in  q,  (qmax),  and  then  a  local  minimum,  (qmm).  We  define  the  depth 


of  the  surface  depression,  d,  as: 
,  max  min 

d  =  q  -q 


(61) 


Figure-14  shows  the  depth  of  the  depression,  d,  as  a  function  of  the  Weber 

number,  We.  There  is  no  depression  (ie,  d=0)  up  to  a  critical  Weber  number 
(Wec).  For  We  larger  than  the  critical  value,  the  depth  of  the  depression  increases 

rapidly  with  the  We. 

The  depression  width  8  is  defined  as  the  gas  gap  corresponding  to  a 
depression  of  d/2.  Figure-15  shows  5  as  a  function  of  We.  We  see  that  for  high 
We,  5  levels  off  to  a  value  of  about  0.013. 

DISPERSION  ANALYSIS 

Figure- 16  shows  a  typical  plunging  liquid  jet.  The  axisymmetric  liquid  jet 
leaves  the  nozzle  with  an  angle  0',  diameter  Dj  and  velocity  ty  The  jet  velocity, 

diameter  and  impact  angle  are  modified  by  gravity;  the  values  at  the  pool  level 
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being  Uy,  D  and  0,  respectively.  The  jet's  inertia  produces  a  meniscus  on  the 

surface  of  the  pool.  It  is  proposed  that  the  gas  entrainment  is  produced  by  an 
instability  of  the  gas/liquid  interfaces.  The  destabilizing  factors  are  the  liquid  jet’s 
velocity  and  gravity. 

The  situation  is  idealized  as  shown  schematically  in  Figure-17.  Because 
the  induced  air  gap  is  very  small  compared  to  the  diameter  of  the  jet  (D),  cartesian 
symmetry  is  assumed.  Region  1  is  the  half  width  of  the  liquid  jet  (hL  =  D/2), 

region-2  is  the  gas  gap,  h2  =  8,  and  region-3  is  the  liquid  pool's  width  (h3  =  <», 
u3  =  0). 

The  following  hypotheses  are  made: 

(i)  The  fluids  are  inviscid,  incompressible  and  the  flow  is  irrotational. 

(ii)  All  interfaces  are  planar. 

(iii)  The  component  of  the  gravity  in  the  y-direction  affects  the  stability  of 

the  flow.  The  effect  of  the  gravity  in  the  x-direction  is  to  accelerate  the 
flow  (i.e.,  to  make  hp  h2  and  hg  functions  of  x).  Because  the  region  of 

interest  in  the  x-direction  is  small,  this  effect  has  been  neglected. 

(iv)  A  sinusoidal  perturbation  of  amplitude  "a"  is  imposed  on  interface 
A.  This  produces  a  perturbation  of  amplitude  "b”  on  interface  B. 
Both  perturbations  travel  with  celerity  C. 

further  assum  that  amplitude  "a"  is  small  enough  in  order  to 
assure  the  validity  of  linear  analysis.  If  we  make  a  Galilean  transformation  with 

velocity  C,  the  perturbations  are  standing  waves  and  the  velocity  of  the  fluids  i  re: 
(ux  -  C),  (U2  -  C),  (u3  -  C),  in  regions  1,  2  and  3,  respectively. 

The  appropriate  conservation  laws  are  given  below: 

Mass  Conservation 

V»ii  =  0  (19) 
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i 


where  we  have  assumed  that  the  fluids  are  incompressible.  Since  the  flow  is  also 
assumed  to  be  invisdd  and  irrotational,  Eq.  (19)  reduces  to: 


V2<j>  =  0 

where,  u  =  -  V<j>,  and  $  is  the  velocity  potential, 
yrnimmtnm  Conservation 


|4+u.Va  =  .ivP+g 


hence, 


H  + 1  1 11 1  2  +  p/p  +  gyy  =  const. 


(20) 


(21) 


(22) 


where, 


u  = 


9x  dy 


dy  dx 


(23) 


Using  eqs.  (22)  and  (23)  and  complex  variables  we  obtain  the  following 
dipersion  relationship: 

[-o12  k2  +  (ux  -  C)2  kpj  coth  (khp  +  (u2  *  C)2  kp2  coth  (kh2)  -  gy(P2  -  Pl>] 
x  [-a23  k2  +  (u3  -  C)2  kp3  coth  (kh3)  +  ^  -  C)2  kp2  coth  (kh2)  -  gy(p3  -  p^] 

=  (u3  -  C)4  k2p2  cosech(kh2)  (24) 

Equation  (24)  gives  the  wave  velocity,  C,  for  a  given  wave  number,  k.  If  the 

solution  of  Eq.  (24)  for  C  is  a  complex  number  with  a  negative  imaginary  part  (ie, 
Cj  <  0)  then  the  system  is  unstable,  since  this  implies  exponential  growth  of  the 

form,  Exp[-Cjkt]. 
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Equation  (24)  is  a  more  general  solution  than  the  problem  of  interest, 
nevertheless  it  can  be  reduced  to  the  plunging  liquid  jet  case  by  noting; 


hj  =  D/2  Uj  =  u 

h3  ~  U3  ■■  0 

°12  =  °23  =  0  P  j/Pl  <<:  * 

Hence,  the  dispersion  relation  becomes: 
(u  -  C)2  k  cothfk  7^  +  g  sin  0  - 


Pl  =  P3  =  Pi 
P2  =  Pg 


=  0 


(25) 


(26) 


The  four  solutions  of  Eq.  (26)  are: 


(27a) 


C 


u  ± 


\ 

-  g  sin0 

) 


1 

k  coth  (~^) 


(27b) 


The  first  two  solutions  in  Eq.  (27a)  are  always  real  and  therefore  they  do  not 
produce  interfacial  instability.  In  contrast,  one  of  the  two  roots  in  Eq.  (27b)  leads 
to  unstable  wave  growth  if, 


k  < 


gp]Sin0 

c 


(28) 


GAS  ENTRAINMENT  RATE 

We  now  have  all  the  elements  necessary  to  calculate  the  volumetric  flow 
rate  of  air,  Qa.  Figure-18  presents  a  schematic  of  the  air  entrainment  process. 

Two  interfacial  waves  of  small  amplitude,  but  wave  number  kd,  grow  as  they 
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move  with  speed  C.  At  a  certain  position  they  touch  each  other  entrapping  a 
volume  of  air  proportional  to  the  shaded  area.  The  shaded  area,  A,  is  given  by: 


o 


The  corresponding  volume  of  entrapped  air  is, 

V  =  icDA  s  jtSDXd 

Finally,  the  volumetric  air  flow  is  given  by: 

Qa  =  ^  nD  XdC  =  kSDC  (29) 


Figure- 19  shows  the  volumetric  flow  rate  of  air,  Qa,  measured  by  McKeogh 
&  Ervine  (1981)  as  a  function  of  the  liquid  jet  velocity,  Uj,  for  a  jet  diameter  of, 

D  =  0.0051  m.  The  jet  turbulence  level  in  these  experiments  was  3%  and  the 
distance  of  the  nozzle  to  the  pool  surface  level  was  0.03  m.  We  note  in  Figure-19 
that  Eq.  (29)  agrees  well  with  the  data.  In  order  to  further  test  the  validity  of  the 
model  we  need  to  know  the  amplitude  of  the  perturbation  (i.e.,  the  apparent 
roughness  of  the  liquid  jet). 

PART-TV:  NUMERICAL  RESULTS 

THE  EULERIAN  LAGRANGIAN  MODEL 

The  momentum  balance  for  the  bubbles  in  component  form  for  a  two-dimensional 
steady-state  flow  can  be  written  as: 

D  y 

Xt.u  (30a) 
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(30b) 


D  * 

-  X  JL  s  u 
Dt  P 


PgVb “dT~  *  P^b(!1/  *  +  Cvmp/Vb  ^  *  V)H<  -  “dT 


ip/AbcD'a/  -  Hgl(^  -ug)+(p^-  pg)Vbg 


+CLP/Vbur  x(Vxu^) 


(30c) 


The  so-called  "dirty  water"  Wallis  drag  correlation  was  used  to  compute  the 
drag  on  bubbles: 

Cd-^  <3» 

Eqs.  (30)  form  a  system  of  Ordinary  Differential  Equiations  (ODEs)  with  the 
position  and  velocity  of  the  bubbles  as  dependent  variables,  and  time  as  the 
independent  variable.  In  this  form  it  is  a  initial  value  problem,  with  the  position 
and  velocity  of  the  bubble  at  time  t=0  the  initial  conditions.  We  have  used  the  well 
known  Tolmien  solution  for  the  liquid  velocity  field.  This  is  a  single-phase 
calculation  of  a  free  jet  using  Prandtl’s  mixing  length  theory  to  model  the 
turbulence.  The  solid  lines  in  Figure  20  show  the  trajectories  of  three  bubbles  with 
different  initial  conditions  for  a  lift  force  coefficient  of  0.5  (solid  lines)  and  a  lift 
force  of  0.25(dashed  lines). 

Figures-20  and  21  show  the  emergence  time  and  the  maximum  depth  as  a 
function  of  the  jet  velocity.  These  are  two  quantities  of  interest  for  mass  transfer 
calculations.  The  emergence  time  is  the  time  that  the  bubble  is  under  water.  We 
see  in  Figure-21  that  the  general  trend  is  that  the  emergence  time  increases  with 
the  jet  velocity.  The  maximum  depth  also  increases  with  the  jet  velocity  up  to  a 
V=2.5m/s.  For  velocities  higher  than  2.5  m/s  the  maximum  depth  is  not  very 
sensitive  to  jet  velocity. 
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THE.TffiQrELUID  MODEL 


Different  phenomenologically-based  models  for  two-phase  flows  have  been 
proposed  in  the  past.  The  drawback  of  many  of  these  models  is  that  they  are  only 
applicable  to  particular  problems.  On  the  other  hand,  mechanistically-based 
models,  commonly  known  as  Two-Fluid  Models  (TFM),  have  been  developed 
[Ishii,  1975;  Delhaye,  1968  and  Drew  &  Lahey,  1979]. 

For  an  adiabatic  plunging  liquid  jet  entraining  air  bubbles,  we  have  the 
following  local,  instantaneous  conservation  equations: 


Mass 


aPk 

-^  +  V*(pkyk)  =  0 


(k  =  g,  1)  (32) 


Momentum 

it  (Pk  +  v#(Pk  *k£k>  =  v*2k  +  Pkl  < k  =  <33) 


where, 

5k  =  '  Pkl +  Jk  ^ 

and,  pg,  pj,  Yg,  vj,  pg,  pj,  Tg,  Xj  are  the  phasic  densities,  velocities,  pressures  and 

shear  stresses  of  the  gas  and  liquid  phases,  respectively. 

These  local,  instantaneous  conservation  equations  may  be  appropriately 
averaged  to  obtain  the  two-fluid  model.  Ishii  [1975]  and  Delhaye  [1968]  have 
proposed  time  and  spatial  averages.  However,  the  ensemble  average  [Drew  & 
Lahey,  1979]  seems  to  be  the  most  fundamental  type  of  averaging.  An  ensemble 
average  is  defined  as, 
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(35) 


f(x,t)=jf(S,t;5)P(5)d5 

s 

were  f  is  the  function  we  want  to  average,  x  is  the  spatial  coordinate,  t  is  time,  £  is 
a  parameter  that  determines  a  particular  realization,  P(£)  is  the  probability 
density  function,  and  S  is  the  set  of  all  realizations.  We  note  that  P(^)  satisfies, 

JP(5)d$«L0  (36) 

s 

Let  us  define  the  phase  indicator  function,  Xk^.t),  such  that  it  is  unity  if  phase-k 
is  present  at  x  and  time  t,  and  is  zero  otherwise. 

Thus,  the  volume  fraction  of  phase-k,  ak,  is, 

«k  *  (37) 

The  physical  interpretation  of  ak  is  as  follows;  ak(x,t)  is  the  fraction  of  all 

realizations  in  which  phase-k  is  present  at  location  x  at  time  t. 

We  refer  the  reader  to  Arnold  [1988]  and  Park  [1992]  for  complete  details  on 
the  derivation  of  the  two-fluid  model.  We  present  here  only  the  final  results. 

The  ensemble-averaged  continuity  equation  for  phase-k  is: 


Mass 


3(akPk) 


dt 


+  v,(°kPkvk)  =  0 


(k  =  g  or  1)  (38) 


where, 


Pk 


xkpk 

ak 
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(39) 


is  the  ensemble-averaged  phasic  density. 


Mgi  - -M«  +  V  •  [ag[os  +  (Pgi  -  )) 


We  note  that,  for  monodispersed  spherical  bubbles,  Laplace's  equation  yields, 

P~-P^  =  2c/Rb  (47) 

where,  a  is  the  surface  tension  and  Rb  the  mean  bubble  radius.  The  bubble's 
surface  stress  tensor,  is  given  by  [Park,  1992], 


2s=p<K^rXr  +  bs(vr*»r)l] 


where,  vr  =  vg  -  \lt  is  the  average  relative  velocity,  and,  using  inviscid  flow 
theory,  the  coefficients  a  and  b  in  Eq.  (48)  can  be  analytically  computed  for 


spherical  bubbles  to  be: 


as  =*20 


h  A 
bs  ~  20 


TURBULENCE  MODELING 

The  total  Reynolds  stress  tensor  for  the  continuous  liquid  phase  is  given  by 
superposition  as, 

Re  Re  Re  . 

il  -Jl(BI)  +il(SI)  (50' 

where,  for  bubbly  two-phase  flows  the  bubble-induced  shear  stress  is  given  by 
[Arnold  1988], 


allKBI)  =  agp[al  ^r— r  +  bl<2r*2r>l] 


where  the  coefficients  aj  and  ty  can  be  analytically  computed  for  spherical  bubbles 
to  be, 
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8l{  =  -  1/20 


b|  =  -  3/20 


(52) 


We  note  that  Iksd  is  the  shear-induced  Reynolds  stress  which  comes  from  the 
classical  k-e  model  [Rodi,  1984]: 


DlCaiW  [  1 

- Dt - *  V lalul  Vkl(SI)J  +  al(PI  •  £1> 


(53a) 


D][oti£l] 


Dt 


=  V* 


Ve, 


+  a^Pg  -  e£) 


(53b) 


where, 


1  _2 
kl(BI)  =  2  ttg  Cvm  1  — r  1 

(54a) 

kl  =  kl(SI)  +  kl(BI) 

(54b) 

,2 

t  1(SI) 

U1  ~  S  ej 

(54c) 

P]  =  Uj  (Vvj  +Vj  V):Vvj 

(54d) 

£jP  1 

pe  =  Cie  k 

kl(SD 

(54e) 

ee  =  C2e  el/kl(SI) 

(54f) 

PrJ  =  1.3 

(54g) 

[Rodi,  1984],  =  0.5478,  Cd  =  0.1643,  Cle  =  1.44,  C2e  = 

1.92  are  the  single- 

phase  flow  k-e  model  coefficients. 

Using  these  results,  the  Reynolds  stress  tensor  is  given  by: 

JKSI)  =  ^1  [v*l  +  ilv]  (55) 
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After  averaging,  we  have  more  unknowns  than  equations.  Thus,  we  need 
to  constitute  some  of  the  averaged  quantities  in  terms  of  the  state  variables,  ak,  vk 

and  pk.  This  closure  process  is  necessary  because  we  have  lost  information  due 

to  averaging  and  we  must  reintroduce  the  essential  physics  which  was  lost. 

Cell  model  averaging  was  successfully  applied  to  a  dilute  mixture  of  liquid 
and  gas  bubbles  by  Arnold  [1988]  and  Park  [1992].  We  divide  the  flow  field  into 
"cells",  each  of  which  have  only  one  spherical  bubble  inside.  Using  inviscid  flow 
theory  we  may  compute  the  pressure  distribution  around  the  bubbles  and  thus 
deduce  the  various  interfacial  forces.  The  assumptions  are  that  both  phases  are 
inviscid,  incompressible  and  have  constant  thermophysical  properties,  the 
bubbles  are  spherical  and  can  be  treated  as  a  dilute  dispersion  of  spheres,  the  non¬ 
uniformities  in  the  distribution  of  the  dispersed  phase  are  small  and  the  velocity 
gradients  of  both  phases  are  small. 

Using  the  results  of  Arnold  [1988]  and  Park  [1992]  the  resultant  two-fluid  phasic 
momentum  equations  are: 


Momentum  Conservation  -  Gas  Phase 

a 


s  (<Wg)  +  v  •  (<w*l*)  =  •  V^i 


^vm^gPl 


Lii+^*V 


CrotagPft-  X  V  x  yg  -  CLpjagvr  x  V  x  ^ 


-  (Cx  +  C2  -  2Cp  -  2bs)  ccgpjyj.  •  VvJ  +  (ag  -  C2)  agp]Vr  •  Vyr 
-  -  CD  _  _ 

+  (as  -  C2)agpj(V  •  yr)  yr  +  agPgg  -  -§-  pj  A"'yr  I  yr  I  -  GroP^Vcig 


(56a) 
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a  _  _  ..  _  2 

ft  (ajp^)  +  V  •  (cqppq  Yt)  =  -  ajVpj  +  (Cp  +  bs  +  bp  pi  I yr  I  Vag 

+  CvmagPl  [(I  +  vg  •  v)  vg  -  +  7,  •  v)  v,  ] 

+  CrotagP^ r  x  V  x  2g  +  CLPlag2r  x  V  x  Yj  +  (C2  +  a1)agp1  yr  (V  •  yr) 

+  (Cx  +  C2  +  2bp  agpiir  •  VyJ  +  (C2  +  ap  agpj  yr»  Vyr 

—  -  C  j-)  _  _ 

+  (as  +  ap  p!  (vr  •  Vag)  yr  +  g.  +  ~g“  Pi  A"' yr  I  Yr  I  +  GrDpjkjVoCg  (56b) 

where  both  phases  were  assumed  to  be  incompressible  and  the  following  values  of 
the  closure  parameters  have  been  used  [Park,  1992]: 


=  - 1/20 

as  =  -  9/20 

bj  =  -  3/20 

bs  =  3/20 

Cp  =  1/4 

Cx  =  5/4 

C2  =  -  9/20 

Gj<d  =  0*1- 

CL=  0.05 

Crot =  005 

For  numerical  purposes  it  is  important  to  rewrite  Eq.  (56a)  as  follows.  For 
the  air/water  flows  under  consideration  Eq.  (38)  can  be  employed  to  show  that  the 
left  hand  side  of  the  gas  phase  momentum  equation  can  be  rewritten  in 
Lagrangian  form  as: 

I  (“gPgV  +  v  •  <<W«V  =  P6°g  m  (S7) 

Grouping  the  right  hand  side  of  Eq.  (57)  with  one  part  of  the  virtual  mass  force  in 
Eq.  (56a)  we  obtain: 


Dgv  DgV 

Pg“g^t  +  CvmPl  “g-Dt 

Notice  that  we  have  kept  the  other  part  of  the  virtual  mass  force  associated  with 
liquid  phase  acceleration  in  the  second  term  on  the  right  hand  side  of  Eq.  (56a).  It 
has  been  found  that  writing  the  gas  phase's  acceleration  as  in  Eq.  (58)  helps 
numerical  convergence. 

Equations  (56)  were  numerically  evaluated  using  the  finite  difference 
formulation  of  Patankar  [1980]  in  the  well-known  PHOENICS  code.  The 
dependent  variables  were  calculated  and  stored  at  discrete  points  on  the  grid.  To 
prevent  pressure  "checker  boarding"  [Patankar,  1980]  a  staggered  grid  was  used. 
The  velocities  are  calculated  at  the  locations.  The  cell  surrounding  point  P  is 
often  called  the  continuity  cell.  The  velocities,  vy  and  vz,  were  computed  at  the 
cell  boundaries  and  the  pressure  and  void  fraction,  p  j  and  ag,  were  computed  at 
the  continuity  cell  center  (P). 

According  to  the  differencing  procedure,  the  conservation  equations  were 
first  integrated  over  the  control  volume  that  surrounds  the  node.  The  resulting 
integrals  were  then  approximated  using  the  nodal  values  and  algebraic  difference 
equations  were  obtained,  where  the  discrete  equations  had  an  implicit 
formulation. 

For  the  boundary  conditions  we  used  the  average  axial  liquid  and  gas 
velocities,  vJz  and  vg2,  the  liquid  turbulent  kinetic  energy,  v^,  and  the  gas  void 
faction  ag  that  were  actually  measured  at  the  inlet  of  the  integration  domain. 
Figure-23a  shows  the  average  axial  liquid  velocity,  vlz.  The  open  circles  are  the 
experimental  values  at  the  inlet  of  the  domain.  The  solid  curve  is  the  computed 
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Vfc  at  the  axial  position  z  =  2.5  mm  (i.e.,  at  the  first  velocity  node),  and  y  =  0.1  m 

corresponds  to  the  planar  jet's  centerplane. 

Figure-23b  shows  the  computed  vlz  at  z  =  31  mm.  The  open  circles  are 

experimental  points.  We  can  see  that  the  agreement  is  quite  good.  The  spreading 
of  the  jet  is  well  predicted  and  the  underprediction  at  the  center  line  velocity  is 
similar  to  that  observed  in  single-phase  flows  [Rodi,  1984].  Figure-23c  shows  the 
computed  vlz  at  z  =  59  mm.  The  open  circles  are  again  experimental  points,  and 

the  trend  is  similar  to  Figure-23b. 

Figure-24a  shows  the  gas  volume  fraction  as  a  function  of  the  lateral 
position.  The  open  circles  are  the  experimental  values.  The  solid  curve  is  the 
computed  ag  at  the  axial  position  z  =  1.25  mm  (i.e.  at  the  first  gas  volume  faction 

node).  Figure-24b  and  24c  show  the  gas  volume  fraction  as  a  function  of  the 
lateral  position  at  distances  from  the  integration  domain  inlet  of  z  =  31  mm  and  z 
=  89  mm,  respectively.  We  can  see  that  the  agreement  is  good,  however,  it  can  be 
noticed  that  the  model  tends  to  overpredict  gas  dispersion.  For  example  in 
Figure-24c  the  experimental  peaks  are  higher  than  the  predicted  ones  and  the 
experimental  center  plane  valley  is  somewhat  deeper  than  predicted. 

PART  -V:  SUMMARY  AND  CONCLUSIONS 

Detailed  measurements  were  taken  of  the  three-dimensional  void  fraction 
and  liquid  velocity  fields  beneath  a  plunging  liquid  jet.  In  particular,  detailed 
new  LDA/FPBA  data  have  been  taken  of  the  air  carryunder  process  associated 
with  a  plunging  cylindrical  liquid  jet.  In  addition,  similar  data  (not  reported 
herein)  was  taken  with  the  planar  nozzle  shown  in  Fig-25.  The  size  distribution  of 
the  entrained  air  bubbles  has  been  measured  directly  with  the  FPDA  system.  It 
was  found  that  for  a  smooth  jet  the  entrained  bubbles  were  very  small  (db  =  120 
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jim).  For  this  case  the  slip  ratio  was  nearly  unity,  and  the  turbulent  intensity  of 
the  liquid  phase  was  comparable  to  single-phase  liquid  jets. 

The  phase  doppler  anemometer  system  was  found  to  be  especially  well 
suited  for  small  bubble  size  measurements  in  two-phase  jets.  With  our  particular 
setup  we  could  measure  bubble  diameters  up  to  about  1  mm.  Because  of  the 
larger  bubbles  present,  we  did  not  use  the  FPDA  to  measure  the  bubble  diameters 
with  the  rough  jet. 

Indeed,  for  a  turbulent  (ie,  rough)  liquid  jet  the  entrained  bubble  sizes  were 
of  the  order  of  d  =  2  mm,  and  the  slip  ratio  was  close  to  calculated  values  based  on 
the  terminal  rise  velocity  of  a  single  bubble.  Moreover,  the  turbulence  intensity  of 
the  liquid  jet  had  two  components,  one  due  to  shear-induced  turbulence  and  the 
other  due  to  bubble-induced  pseudo-turbulence.  Both  components  of  turbulence 
were  of  the  same  order  of  magnitude. 

The  FPDA  system  was  found  to  be  especially  well  suited  for  small  bubble 
size  measurements  in  two-phase  jet  flows. 

A  state-of-the-art  two-fluid  model  obtained  using  ensemble  averaging  has 
been  derived  and  was  closed  using  cell  average  model.  This  approach  provides 
equations  for  multiphase  flows  that  are  mechanistically-based  (as  opposite  to 
empirical).  The  rigorous  derivation  of  the  cell  average  model  provides  exact 
constitutive  equations  for  the  inviscid  limit.  One  does  not  expect  the  values  of  the 
constants  from  cell  model  averaging  to  be  correct  for  very  viscous  flows  but  they 

provide  a  good  framework  from  which  to  start.  In  particular,  it  is  known  [Wang 
et  al,  1987]  that  the  lift  coefficient,  CL,  decreases  as  liquid  viscosity  increases.  In 

this  study,  a  lift  force  coefficient  of  CL  =  .05  has  been  used  instead  of  the  theoretical 
inviscid  limit  value  of  CL  =  .25.  All  other  parameter  values  used  in  this  work 
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corresponded  to  the  inviscid  limit  values.  The  agreement  with  the  experimental 
data  was  remarkable. 

The  k-s  model  seems  to  be  adequate  for  this  calculation,  however  the 
observed  differences  in  phasic  dispersion  indicate  some  inadequacies  in  the 
turbulence  modeling  which  should  be  considered  in  future  studies. 

This  report  also  presents  the  results  of  an  analysis  based  on  a  non-singular 
perturbation  technique  for  the  surface  depression  produced  by  a  plunging  liquid 
jet.  The  first  three  terms  of  a  Taylor  series  expansion  in  the  Weber  number  have 
been  used  to  give  the  pool  surface  depression.  This  approximation  of  the  surface 
depression  gives  correct  values  for  small  and  moderate  Weber  numbers  (ie,  We  < 
1).  However,  for  We  greater  than  unity,  higher  order  terms  may  become 
important  and  the  analysis  presented  herein  is  no  longer  valid.  The  Bond 
number  appears  as  a  parameter  in  the  system  of  equations,  thus  no  expansion  is 
necessary  for  the  Bond  number.  Indeed  the  analysis  presented  herein  is  valid  for 
all  Bo. 

The  results  described  in  this  report  show  that  as  the  Weber  number  is 
increased,  the  terms  that  gain  importance  (i.e.,  the  terms  of  higher  order) 
correspond  to  a  surface  depression  that  is  increasingly  narrower  in  the  horizontal 
direction.  For  a  Weber  number  of  the  order  of  unity  the  surface  tension  is  no 
longer  able  to  keep  the  system  stable  and  an  Helmholtz  instability  leads  to  air 
entrainment  [Bonetto  et  al,  1993].  That  is,  for  values  of  the  surface  tension  going 
to  infinity  (i.e.,  Weber  number  going  to  zero),  the  slope  of  the  surface  is  very  small, 
however,  as  Weber  number  is  increased  the  slope  also  increases.  For  a  critical 
value  of  the  Weber  number,  the  slope  of  the  surface  is  such  that  the  surface 
tension  is  not  able  to  keep  the  deformed  pool  surface  away  from  the  plunging 
liquid  jet  and  air  entrainment  is  produced. 
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We  have  assumed  in  this  report  that  the  position  of  the  interface  can  be 
written  as  q(x)  ,  eq.  13a,  that  is  a  steady  state  exists  where  the  position  of  the 
interface  is  a  single-valued  function  of  the  lateral  position.  For  the  critical  value 
of  Weber  number  that  corresponds  to  the  threshold  for  air  entrainment  the 
interface  is  not  single-valued.  Moreover,  it  is  likely  that  the  air  entrainment 
process  involves  transient  phenomenon.  Thus  we  believe  that  our  analysis  is  not 
really  appropriate  to  compute  the  threshold  value  of  Weber  number  for  air 
entrainment.  However,  we  are  able  to  describe  the  route  to  air  entrainment 
qualitatively. 

We  have  seen  that  the  higher  the  order  of  the  position  of  the  interface 

iterate,  the  larger  the  maximum  slope  of  the  iterate.  For  small  Weber  numbers 

the  low  order  iterates  dominate.  As  Weber  number  is  increased,  the  slope  of  the 

interface  is  increased  as  the  higher  order  terms  dominate.  For  some  threshold 
value  of  Weber  number,  Wet,  the  interface  has  a  very  large  slope  (ie,  the  tangent  of 

the  interface  is  almost  vertical).  For  We>  Wefc  the  inertia  forces  overpower  the 

capillary  forces  and  the  interface  is  no  longer  stable. 

The  results  presented  in  this  report  cannot  predict  all  details  of  the  route  to 
the  pool  surface  instability  that  produces  air  entrainment  for  two  reasons.  First, 
it  is  not  clear  that  the  approximation  for  the  surface  depression  is  valid  for  Weber 
numbers  that  are  high  enough  to  produce  air  entrainment.  Second,  the  shape  of 
the  surface,  T](x),  has  to  be  a  single-valued  function  of  x,  however  it  appears  that 
at  the  point  where  air  entrainment  commences,  for  every  x  (horizontal  position) 
there  is  more  than  one  rj  (vertical  position  of  the  surface  depression). 

It  appears  that  it  would  be  useful  to  compute  the  surface  position  using  an 
appropriate  multidimensional  Computational  Fluid  Dynamics  (CFD)  tool  having 
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surface  tracking  capability.  This  will  allow  the  relaxation  of  the  restriction  on 
Weber  number. 

This  report  presents  the  results  of  a  linear  stability  analysis  for  an  invisdd, 
irrotational  and  incompressible  liquid  jet  impacting  a  static  pool  of  liquid.  It  has 
been  shown  that  a  Helmholtz-Taylor  instability  can  occur  and  that  this  instability 
may  lead  to  air  entrapment  and  carryunder.  Moreover,  it  is  shown  that  the 
thickness  of  the  annular  air  gap  is  not  a  function  of  liquid  jet  velocity,  however 
increasing  the  turbulence  level  of  the  liquid  jet  increases  the  wave  action  on  the 
surface  of  the  jet  which  leads  to  increased  air  entrainment. 

It  appears  that  it  would  be  useful  to  have  future  experimental  research  on 
this  topic  focused  on  the  effect  of  liquid  jet  viscosity,  nonsymmetric  situations,  and 
nonlinear  phenomena.  Moreover,  it  should  be  useful  to  couple  this  type  of 
analysis  with  the  induced  two-phase  jet  spreading  problem  below  the  surface  of 
the  pool.  This  is  probably  best  done  using  appropriate  interface  tracking 
computational  fluid  dynamic  (CFD)  tools. 
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Figure  1  Schematic  of  the  Experimental  Apparatus 


Figure  2  Schematic  of  the  Conical  Nozzle 
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Figure  4  Liquid  Velocity  Histogram  -  Rough  Jet 

(z  =  33  mm;  w*  =  0. 144  kg/s;  h  =  29.9  mm;  r  =  0) 
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Figurc-9:Void  fraction  as  a  function  of  the  lateral  position 
for  different  axial  positions. 


Figure  11  Illustration  of  the  surface  depression  produced  by  a  plunging 
liquid  jet. 


Figure  12  Iterates  of  the  surface  position,  Tij(x),  corresponding  to  Bo-0 
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Figure  13  Shape  of 
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Fig.  16  The  Plunging  Liquid  Jet. 


* 

•.  ..  -  1  . 

r  ss  u  ‘ 
•—  f  •»  i  r--3 


Fig.  17  Schematic  Diagram  of  the  Plunging  Liquid  Jet  Just 
Below  the  Pool  Surface. 
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Fig.  18  Mechanism  Responsible  for  Air  Entrainment. 


CL  =  0.25 
CVM  =  0.5 
6=45 
V=5m/s 
Rb=0.002m 


Figure  20  Bubble  trajectories  for  Cl=0.25  and  different  initial  positions 

a)  y(t=0)=-h/2 

b)  y(t=0)=0 

c)  y(t=0)=+h/2 

The  dashed  lines  are  the  base  solutions 
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Figure  21  Emergence  time  t  as  a  function  of  the  jet  velocity  for  different 

initial  positions: 

a)  y(tas0)=!-h/2 

b) y(t=0)=0 

c)  y(t=0)=+h/2 


dm(m) 


Figure  22  Maximum  depth  dm  as  a  function  of  the  jet  velocity  for  different 
initial  positions: 

a)  y(t=0)=-h/2 

b)  y(t=0)=0 

c)  y(t=0)=+h/2  . 


Axinl  velocity  ns  a  function  of  the  lateral  coordinate  (y  =  0.1 
corresponds  to  the  jet  center  plane)  -  Inlet  conditions. 


Void  fraction  ns  n  function  of  the  lateral  coordinate,  (y  «  0.1 
corresponds  to  the  jet  center  plane.  The  circles,  o,  are  experimental 

results) 


25  Schematics  of  the  test  section. 


